function [xkesi,ykesi,xyita,yyita,J] = xykesiyita(x,y,dkesi,dyita,nelx,nely)
%XYKESIYITA Summary of this function goes here
%   Detailed explanation goes here
xkesi = zeros(size(x));
ykesi = zeros(size(x));
xyita = zeros(size(x));
yyita = zeros(size(x)); % solve for x_yita, x_kesi, y_yita, y_kesi

for i = 2 : nelx-1
    for j = 2 : nely-1
    xkesi(i,j) = (x(i+1,j)-x(i-1,j))/(2*dkesi);
    xyita(i,j) = (x(i,j+1)-x(i,j-1))/(2*dyita);
    ykesi(i,j) = (y(i+1,j)-y(i-1,j))/(2*dkesi);
    yyita(i,j) = (y(i,j+1)-y(i,j-1))/(2*dyita);
    end
    xkesi(i,1) = (x(i+1,1)-x(i-1,1))/(2*dkesi);
    xyita(i,1) = (x(i,2)-x(i,1))/dyita;
    ykesi(i,1) = (y(i+1,1)-y(i-1,1))/(2*dkesi);
    yyita(i,1) = (y(i,2)-y(i,1))/dyita;
    
    xkesi(i,nely) = (x(i+1,nely)-x(i-1,nely))/(2*dkesi);
    xyita(i,nely) = (x(i,nely)-x(i,nely-1))/dyita;
    ykesi(i,nely) = (y(i+1,nely)-y(i-1,nely))/(2*dkesi);
    yyita(i,nely) = (y(i,nely)-y(i,nely-1))/dyita;
end

for j = 2 : nely-1
    xkesi(1,j) = (x(2,j)-x(1,j))/dkesi;
    xyita(1,j) = (x(1,j+1)-x(1,j-1))/(2*dyita);
    ykesi(1,j) = (y(2,j)-y(1,j))/dkesi;
    yyita(1,j) = (y(1,j+1)-y(1,j-1))/(2*dyita);
    
    xkesi(nelx,j) = (x(nelx,j)-x(nelx-1,j))/dkesi;
    xyita(nelx,j) = (x(nelx,j+1)-x(nelx,j-1))/(2*dyita);
    ykesi(nelx,j) = (y(nelx,j)-y(nelx-1,j))/dkesi;
    yyita(nelx,j) = (y(nelx,j+1)-y(nelx,j-1))/(2*dyita);
end
    xkesi(1,1) = (x(2,1)-x(1,1))/dkesi;
    xyita(1,1) = (x(1,2)-x(1,1))/dyita;
    ykesi(1,1) = (y(2,1)-y(1,1))/dkesi;
    yyita(1,1) = (y(1,2)-y(1,1))/dyita;
    
    xkesi(1,nely) = (x(2,nely)-x(1,nely))/dkesi;
    xyita(1,nely) = (x(1,nely)-x(1,nely-1))/dyita;
    ykesi(1,nely) = (y(2,nely)-y(1,nely))/dkesi;
    yyita(1,nely) = (y(1,nely)-y(1,nely-1))/dyita;
    
    xkesi(nelx,1) = (x(nelx,1)-x(nelx-1,1))/dkesi;
    xyita(nelx,1) = (x(nelx,2)-x(nelx,1))/dyita;
    ykesi(nelx,1) = (y(nelx,2)-y(nelx,1))/dkesi;
    yyita(nelx,1) = (y(nelx,2)-y(nelx,1))/dyita;
    
    xkesi(nelx,nely) = (x(nelx,nely)-x(nelx-1,nely))/dkesi;
    xyita(nelx,nely) = (x(nelx,nely)-x(nelx,nely-1))/dyita;
    ykesi(nelx,nely) = (y(nelx,nely)-y(nelx-1,nely))/dkesi;
    yyita(nelx,nely) = (y(nelx,nely)-y(nelx,nely-1))/dyita;
 
    J = xkesi.*yyita - xyita.*ykesi;   %approximate for J


end

